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Chapter 1 
Introduction 


A common method of analysis for high-frequency electromagnetic and acoustic scat- 
tering and diffraction problems involves the use of radiation integrals as well as plane 
wave integral representations for the fields, with the asymptotic approximations to 
the various scattering mechanisms found from the critical point contributions of 
the integrand. The incomplete Airy integrals [1] serve as canonical functions for 
the uniform asymptotic approximation of a class of integrals characterized by two 
stationary phase points that are arbitrarily close to one another or to an integra- 
tion endpoint [2]. Integrals of such analytical properties describe transition region 
phenomena associated with composite shadow boundaries resulting from the conflu- 
ence of two stationary points and an endpoint in the integration interval. A typical 
problem of particular interest, where transition region phenomena of this type exist, 
involves the scattering from smoothly indented boundaries containing an edge as 
illustrated in Figures 1.1 and 1.2, When the reflection shadow boundary (RSB) is 
not in the immediate vicinity of the smooth caustic, the conventional UTD edge 
diffraction coefficient [3] which involves the Fresnel integral as a canonical function 
can be used to effectively describe the edge diffracted field behavior in the neighbor- 
hood of the reflection shadow boundary. Furthermore, the ordinary Airy integrals 
and their derivatives are the appropriate canonical functions for the description of 
Uniform Geometrical Theory of Diffraction 
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Reflection Shadow Boundary (RSB) 

Figure 1.1: Scattering and diffraction from a concave boundary containing an edge. 

the high-frequency fields in the neighborhood of the smooth caustic [4]. However, 
when there is a confluence of both reflected and caustic type shadow boundaries 
as shown in Figures 1.1 and 1.2, neither the Fresnel integral nor the ordinary Airy 
integral adequately describe the transition region phenomena, and they must be 
appropriately replaced by the incomplete Airy function. 

The formulation of uniform ray optical solutions for problems involving compos- 
ite shadow boundary transition phenomena, that are useful for engineering purposes, 
requires an efficient and accurate method for computing the incomplete Airy func- 
tions. In the original work of Levey and Felsen [1], several other diffraction problems 
whose solutions involve the incomplete Airy functions are also described, and some 
general and asymptotic characteristics of these functions are examined. However, 
the asymptotic formulae given in [1] are valid when the saddle points are sufficiently 
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Reflection Shadow Boundary (RSB) 


Figure 1.2: Scattering and diffraction from a concave-convex boundary containing 
an edge. 


far apart or far removed from the integration endpoint. When the two saddle points 
move close together or coalesce to form a caustic, the asymptotic formulae break 
down and an alternative computation method for the incomplete Airy functions is 
needed. 

In this report, a convergent series solution for the incomplete Airy functions is 
rigorously derived, and asymptotic expansions involving several terms are also devel- 
oped for large argument approximations. The higher order terms in the asymptotic 
expansions greatly improve the accuracy of the asymptotic formulae, and also im- 
prove computational efficiency by limiting the region of the argument space where 
the more time consuming series solution form should be used. Thus, the combination 
of the series solution form with the more accurate asymptotic formulae provides for 
an efficient and accurate computation of the incomplete Airy functions for the entire 
range of their arguments. These results would allow the formulation of uniform ray 
optical solutions for high-frequercy scattering and diffraction problems that are use- 
ful for engineering purposes. A method for uniformly evaluating certain stationary 
phase integrals using the incomplete Airy integral as a canonical function is briefly 
discussed in Appendix A. The details of a systematic uniform asymptotic analysis 
for particular applications that involve transition region phenomena describable by 
the incomplete Airy functions will be reported a separately. 

The outline of this report is as follows: In Chapter 2, some general properties of 
the incomplete Airy functions are reviewed, and a convergent series solution form is 
rigorously derived using the parabolic differential equation. The asymptotic formulae 
for large argument approximations are derived in Chapter 3 using the integral forms 
of the incomplete Airy functions. In Chapter 4, some indicative numerical results 
are presented and discussed, with their accuracy demonstrated via comparison with 
data obtained from direct numerical integration of the integral forms. The regions 
of validity for each formula used in the computations are also shown in this chapter, 
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and some error assessments for the asymptotic results are provided. Finally, the 
main results and accomplishments of this work are summarized in Chapter 5. 
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Chapter 2 

Derivation of the Series Solution 
Form 


In this chapter, a convergent series solution form for the incomplete Airy functions 
is derived. Before proceeding with the derivation, however, some of their general 
properties are briefly reviewed. 

The incomplete Airy functions are functions of two variables and they satisfy the 
parabolic partial differential equation applied by Fock [5] to the study of fields near 
the surface of a smooth convex scattering surface; i.e., 

giW,0 = 0; t = 0,1,2. (2.1) 

In integral form, the solutions of (2.1) are given by: 

g ,(0,O=f eM^Wdz-, * = 0,1,2 * (2.2) 

where the upper limit lies within one of the three sectors of the complex z-plane in 
which the integral converges; i.e., 

2 £ < V’t < (2* + l)^ ; « = 0,1,2. (2.3) 

The contours of integration for the incomplete Airy functions are shown in Figure 2.1. 

In this report, we only examine go(/3,£) since the other two functions, namely 
gj and g 2 , can be obtained from g 0 and the well known complete Airy functions [6]; 


JL-p-jl 

d/3 2 p J d(\ 
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Figure 2.1: Contours of integration for the incomplete Airy functions. 
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gi (0, 0 = &)(£» 0 - 2jrAi(/?) , 


(2.4) 


i.e., 
and 

gaC/^O = go(0,£) - 7r [Ai($) + J*Bi(/3)] , (2.5) 

where 

Ai (j3) = ^-[ e M z+x3/ Vdz, (2.6) 

27T */ 

and 

Bi (0) = J-[ e ii0z+zS/3) dz. (2.7) 

27t y z,2 +^>3 

The contours of integration for the complete Airy functions are shown in Figure 2*2. 
Also, the quantities (3 and £ will be taken as real since in most practical applications 
real /3 and f are of primary interest. However, this is not a requirement for the 
analysis that follows and the resulting formulae are valid for arbitrary values of /3 
and £. In addition, £ will be restricted to positive values since for negative values of 
g 0 may be obtained using the expression: 

go(^> — l£l) = 2wAi(/?) — &)(/?» |£|) > (2.8) 

with (*) denoting the complex conjugate operation. 

In order to obtain a series solution form for gy(/3,£), we begin with the parabolic 
differential equation and assume two independent solutions of the form: 

= E«.«r. (2-9) 

n=0 

OQ 

and y 2 (/?,0 = £ MO/ 3 " • (2.10) 

n=0 

Substituting (2.9) and (2.10) in (2.1) we obtain the following expressions: 

OO OC OC 

^ n(n - l)a„(()|3 n ' 2 - ^ a n (()i3 n+1 - j ^ o' (()^ n = 0, (2.11) 

n = 2 n= 0 n=0 

OO OO OO 

and £n(»-l s 0. (2.12) 

n=2 n=0 n=U 
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For Equations (2.11) and (2.12) to be satisfied, the sum of coefficients of like powers 
of must be zero for any value of £. To obtain the first independent solution, y j , 
we let a](£) = 0, and setting the sum of coefficients of like powers in (2.11) equal to 
zero we obtain: 

o,(0 = o, V ( , (2.13) 

o 2 (0 = |oi(0, (214) 

and o„(0 = + , n > 3 . (2.15) 

n(n — 1) 

Similarly, for the second independent solution, y 2 , we let 6o(£) = 0 and setting the 
sum of coefficients of like powers in (2.12) equal to zero we obtain: 

M0 = 0, VO (2.16) 

M0 = VO (2.17) 

and M0 = ^-alO + AMO 71 > 3. (2.18) 

n(n — 1) 

Thus, y t and y 2 are given by: 

oc 

y,(/3,0 = o„(0 + E “»(0/3" , (2-19) 

n= 2 

and y 2 (/3,0 = MO0 + E MO0" , (2 20) 

n=2 

where a n (0 and b n (() for n > 2 can be expressed in terms of a o (0> M0 anc * t ^ ie ’ r 
derivatives, respectively, via Equations (2.13)— (2. 18). Now, g 0 (/3,0 must be equal 
to the sum of the two independent solutions; i.e., 

DO 

&,(/J,0 = MO + MO/3 + EMO + MO]0 n , (2.21) 

71 = 2 

Q OO 

and ^g 0 (/3 ? 0 = M0 + X) n [°"(0 + MO]#" -1 • (2.22) 

°P 71=2 

Although the coefficients a„(0 an ^ M0 can be combined into a single coefficient, 
keeping them separate greatly simplifies their evaluation in a computer code. 
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Finally, it remains to find expressions for a () (£) and b\(() in order to complete 
the solution. This can be accomplished by applying the proper boundary conditions 
at /? = 0 using the integral form of g(j(/3,£) given in Equation (2.2); i.e., 

fOC . . 

<■«({) = g»(/3. ORo = / e" ' 3 iz , (2.23) 

and 

(0 = Olteo = ( 2 - 24 ) 

The functions a 0 (£) and &i(£) can be expressed in terms of the incomplete Gamma 
function [7], and their computation is straight forward. Details are provided in 
Appendix B. 

When ( — * — oo, our solution should reduce to the series solution form for the 
complete Airy function, Ai (/?) [6]. In this case we have: 

a 0 ({ — » — oo) = 2irAi(0) = 2.23070703 , (2.25) 

6i(£ — > — oo) = 2wAi , (0) = —1.62621025, (2.26) 

and using Equations (2.13)-(2.20) we obtain: 

g„03, {--<*,) = 2jrAi(0) (l + + -^j-0 6 + ~~$r~0 9 + ' * *) 

+ 2^(0)^ + ^ + ^+^^ + ...) 

= 2,rAiQ3), (2.27) 


which is a necessary condition for the validity of our result. 


Chapter 3 


Derivation of Asymptotic 
Formulae 


In this chapter, a pair of asymptotic formulae for the incomplete Airy function, 
9 o(P,£)i involving several terms are derived and serve as large argument forms. We 
begin by introducing the large parameter fl in the integral form of g 0 (/?,£)> and 
examine the integral: 

/„(« r,75«)= r^^dz. (3.1) 

Our objective is to obtain asymptotic expansions for 7 0 as fi — ♦ oo for various 
dispositions of the saddle points and endpoint. Then, the asymptotic formulae for 
9 u{ 0 it) may be obtained using the expression: 

S „(/3,{) = n l / 3 7 0 (^ = /?n - 2 / 3 ,7 = {!5-'/ 3 ;n). (3.2) 

Since we are interested in real /? and £ with £ > 0, the analysis that follows will be 
restricted to real tr and 7, with 7 > 0. 

3.1 Asymptotic Formula for £ > |/3|2 

This case corresponds to the endpoint being far removed from the possibly neigh- 
boring saddle points of the integrand in Equation (3.1), or 7 |<r|2. Although the 

sign of <T is irrelevant in this case, the saddle points z it 2 = ±(— cr) 1 / 2 are taken as 
real for simplicity. Also, the original integration path Pd in Figure 3.1 is deformed 
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into the steepest descent path, leading away from the endpoint at z - 7 and 

into the sector 0 ^ arg(z) < 7 t/ 3 of the complex z-plane. The asymptotic evaluation 
of Jo is then performed using repeated integration by parts that yields: 

27 




~ e 


,jn(<r7+~) 9 /3) 


fiU 

I $7 [<r -f ' 


+ fi 2 


(a + 7 2 Y 




2 j 

12J 7 2 

1 

407 

1207 3 

L(o- + 7 2 ) 4 

(<x + 7 2 ) 5 J 

+ fi 4 

l(. + V) 6 

(o' + 7 2 ) 7 ] 


+ 


ft 5 


» J ) * l / J L V 1 t / V- ' f / J 


[(<r + 7 2 ) 7 (cr + 7 2 ) 8 (<r + 7 2 

Now, using Equations (3.2) and (3.3), the asymptotic formula for gg(/3,£) when 
£ > |/9 1 1/2 is given by: 


g „(£,«) - 

40£ 


[/? + £’ 'V + ^) 3 ' r (,3 + ( 5 ) J (0 + f ! ) 5 

120£ 3 - 40 j 840 j£ 2 1680j£ 4 




+ 


2j 


12 j 


{P+ey (p + i 2 ) 1 {P+ey (P + ( 2 ) 9 \ ' 

The derivative of g 0 (/?,£) is given by: 


dp 


sMS) * 


-( 


+ 


2jV 


H 


8 j - 12£ 60j + 40j£ 2 280£.+ 120j£ 3 


(/ 3 +n 5 (/3+a e 


(£ + £ 2 ) 7 


1680f - 280 j 6720j*e 2 - 168 0^ _ 15120#' 1 


(/? + £ 2 ) e 


(0 + £ 2 ) 9 


(P + ( 2 ) 
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(3.4) 


[p+( 2 (P+t 2 ) 2 (P+rr (P+c 2 ) 4 


(3.5) 


3.2 Asymptotic Formula for /3 <C — 1 


This case corresponds to real and widely separated saddle points (<x «C — 1) in the 
integrand of Equation (3.1), with the endpoint 7 arbitrarily close to the saddle point 
zi = (-o’) 1 as shown in Figure 3.2. For an asymptotic evaluation of Equation (3.1) 
that holds uniformly as the endpoint 7 approaches the saddle point Z\ , we make the 
following transformation [1]: 

q(z) = crz + z 3 /3 = q(zi) + s 1 = — -(— <r) 3 / 2 + s 2 = r(s) , (3.6) 

13 



14 


with arg(.s) restricted so that I, ] ((r,'y,Q) converges as 3 — ♦ 00 . Hence, employing 
(3.6) in Equation (3.1) we have: 


7o(<r, 7 ; n) = r /Me*’ 1 ’ d ’ , 


(3.7) 


with the upper limit taken in the sector 0 < arg(s) < 7r/2 of the complex s-plane. 
The quantities £ and f(s) are given by: 


c = ± 


(T7 + 7 3 /3 + -(-<t) 3/2 


lV2 


; 7 <(-*-) 1/2 > 


and 


f ( Q \ _ ii _ Tlfl - 2s 

f\ 3 ' Jo ~K-,\ ~ 


ds g'(z) <r + z 2 


(3.8) 


(3.9) 


Equation (3.7) can be written as follows [8]: 


n> = e-m-^ {/(«) j~ + 


1 f°° \f(s)~ f(0) 


i-e*"’ d, 

as 


}■ 

(3.10) 


and using integration by parts in the second integral we get: 


| tt™°) - c 


no - m 




where 


F 9< - 3)eia "' ds } ’ 

»/'(«) - /(») + /(o) 


9 < s ) = 


(3.11) 


(3.12) 


In a similar way, the second integral in Equation (3.11) can be written as follows: 


j s(»)e ,n '’ is = g(O) ~j 




2 jil 


9(0 -9(0) 


C 


jo? 




2;Q 


where 


M s ) = 


V( a ) - s(*) + ^(0) 


(3.13) 


(3.14) 
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I==m 


Thus, using (3.13), Equation (3.11) becomes: 


/o(«r,7;0) = f{0) - 2 ^ 0 ( 0 )] ds 


_j_ e >n(<T7+7 3 /3) | 

r i 

[ 2 jn 


l 

l(2jfl )' 


m - m] 


c 


1 

ml h(,) 


+ 


ds 


m ) 2 [ c 


<7(0 -s(0) 


(3.15) 


The same procedure is repeated once more for the last integral in Equation (3.15); 
i.e.. 


/ h{*y 


^ ds 


hm ~m HZ' 


e }U>7 ds 


+e 


>nc : 


l 2jn 


MO - *( 0 ) 


where 


c 


*(*) = 


+ 


1 [ fe(0 _ fe(0) j + 0(n _ 3) ) ( 3 16) 


(2jfi) 


sh'(s) — h(s) + h( 0 ) 


(3.17) 


and finally combining Equations (3.15) and (3.16), the asymptotic expansion of 
Jo(<r, 7 ;f 2 ) when cr «C — 1 is given by: 


Ju( 0 -, 7 ;n) ~ e <t)3/2 

_j_ e in(ff-H-7*/3) 


«°> - zjn* <•> + - (yW*<°>] jf e ’ n " * 

f l r/(o-/(o)i 

1 2jf 


(2jH)3 [ C 


2jf2 

MO - *( 0 ) 


1 [ g(Ojvg(0) ' 

( 2 jH ) 2 l C . 


+ 


w 


^ | + 0(0~ 5 ) . (3.18) 


In order to express the functions /, g , h, and k in Equation (3.18) in terms of <r 
and 7 , we need to derive an expression for f(s) and its derivatives when s & 0. This 
is done using a procedure introduced by Erdelyi [9] that yields: 


/(-) 


and 


d h 

ds k 


/(•) 


1 r(3n/2 - 1)(-1)-' , 

^ (n - l)!r(n/2)(3(-»)W]— ’ 

1 ” T(3n/2 - 1)(-1)""' 

( m T> .S, (« - * - l).r( n /2)[3(— ff) 3 / 4 ]" -1 


(3.19) 

(3.20) 
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Hence, using Equations (3.9), (3.12), (3.14), (3.17), (3.19), and (3.20) we have 

2C_ 

7 2 ’ 

1 


/(C) = 


tr + 7 2 


fn _ ~8tC , 

9K ^ } ~ (o- + 7 2 ) 3 C 2 (-^) ,/4 ’ 

MO = ^4,+ 96,2 


(tr + r-y (<r + 7-) s O(-ff) 1 ' 4 24 0(-ff) 7/4 ’ 


MO = 


/( 0 ) = 


6407 


19207 s 


(tr + 7 2 ) 6 (<r + 7 2 ) 7 + C 6 (-^) 1/4 8C 4 (-o-) 7/4 

77 

+ 3456C 3 (-^) 13/4 ’ 

1 


15 


(-*) 1/4 ’ 
-1 


3 (-«■) ’ 


no) 

5 

2 

24(— cr) 7 / 4 ’ 

no) _ 

-8 

3 

' 27(-<r) 5/2 ’ 

/< 4 >(0) 

77 

8 

" 3456(— tr) 13 / 4 ’ 

/ (5, ( 0) 

-56 

15 

" M(-») 4 ’ 

/**>( 0) 

12155 

48 

“ 82944 (~<r) 19 / 4 

/ (7 >(0) 

-640 

105 

“ 243(-£r) n / 2 * 


and 


( 3 . 21 ) 

( 3 . 22 ) 

(3.23) 


(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 


Also, the Fresnel integral in Equation (3.18) may be expressed in terms of the Fresnel 
transition function F(x ) [3]; i.e., 

n*/ 9 \ 1 

= Q i/2 A(rj), 


/; 


e jna ' ds = fi ,/2 




(3.33) 


where f/(x) is the Heaveside unit step function, and rj = fi 1/,2 C. The general prop- 
erties of the Fresnel transition function and details on its computation are provided 
in Appendix C. Now, using Equations (3.21)— (3.33), the asymptotic expansion of 
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Jo(<7,7;^) when a <C -1 is given by: 

_jln(-<r) 3 / 2 

~ y^z^ s (^)m + ^ wmE ^n,i-,») + 0^- s ), (3.34) 

where 


S{<T' } ft) = 1 — 


Si 8 2 

2 jfl(-<r)3/2 + (2jft) 2 (-o-) 3 


J3 

(2jft) 3 (-<r) 9 / 2 ’ 


(3.35) 


^. 7 ,x;n) = — 


1 

4 

-87 

( 2 jny 

(<T + 7 2 ) 3 

1 

-16 

(2jft) 3 

(<T + 7 2 ) 4 

1 

| 

640 7 

(2j ny 

[(<r + 7 2 ) 6 


2 ft 1 / 2 ‘ 

a + 7 2 ?7( — c r )'^‘ l . 

ft 3 / 2 _ Q 1/2 j», ' 

77 3 ( — cr) 1 / 4 ^(— c) 7 / 4 

967 2 3ft 5/2 n 3/2 j»i _ ft 1/2 a 2 

(a + 7 2 ) 5 J7 S (— o') 1 / 4 ) tj 3 (— o') 7 / 4 »j(— o) 13 / 4 

1920 7 3 15ft 7/2 3ft 5/2 Ji n 3/2 a 2 

(<7 + 7 2 ) 7 + 7J 7 ( — O') 1 / 4 ) “ t? 5 (-o) 7 / 4 + 7/ 3 (-o) 13 / 4 


ft 1/2 * 3 ' 

T/( — 0') 19 / 4 > ’ 


(3.36) 


5, = 0.20833333, s 2 = 0.33420139, and s 3 = 1.02581260. 

The asymptotic formulae for g 0 (/3,£) and its derivative when (3 <C — 1 are then 
obtained using Equations (3.2) and (3.34)— (3.36); i.e., 


goGs.fl - e 7 5 l~'?/r m AM + , 


e-iK-W 


(3.37) 


^ go (/?,0 - JZpyJ* 


{M /2+ iF 


0 ) 


d 


simi) + m-ysw 


+S(0) 




+ e 


jiM+i*/*) 


d 




(3.38) 


where 


s(0) 



2j( rt) , ! T M J (-^) 3 

1.5S! t 3 s 2 

2i(-/3) s/J + (2j) 2 (-iff) J “ 


8 3 

(2j) 3 (— /3) 9 / 2 ’ 


4.5s 3 

(2j) 3 (— ^) n / 2 ’ 


(3.39) 

(3.40) 
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e(P,(,v) = - 


d 

dp 


E(P,t,v) 


dp 


and 


2j [P + t 2 v{ 

1 

•Uj* 

00 

i 

(2j) 2 

[(P + t 2 Y 

1 

-16 

(2j) 3 

i(/3 + £ 2 ) 4 

1 

Ms 

0 

<X> 

1 

(2;) 4 

l(P + ( 2 Y 


+ 


+ 


96£ 2 3 


*1 


*2 


( 2 jy i(p +ey T (p + ey v 5 (-py /4 v 3 (-py /4 v(-pv 3/4 \ 

1920£ 3 , 15 3aj , s 2 


+ 


+ 


S3 


n(-Py 9/4 \ ’ 


(3.41) 


2 j 


-2 


1 


+ 


* - (- py > 2 


XP + t 2 ) 2 H~PY /4 2v 3 (-Py 14 J 

, 1 f 24£ 1 3[4 - C-/g) 1/2 ] 75 t 

(2 jf ((/? + £ 2 )< + 4 tP{-PY' 4 2 t f{-py/* ir,(-P)"/4 

5i[e-(-/3) 1/2 ]) 1 f 64 480f 3 

2 7? 3(- / 3)7/4 / ( 2 jy\(p + eY (P + ( 2 Y *v 5 (-PY /4 

, 15K - (~/?) 1/2 ] .7-1 ~ (-/g) I/a T 13^2 

"* n ** / /^\i / < ( 


Mv) = -ts< — 


2i, 7 (-/3)'/ 4 V(-^)"/ J 2r,-‘(-0y/‘ 4 1 (-0)< 7 I‘ 

» 2 lf-(~^) 1/2 ) ) . 1 3840 ; 13440{ 3 15 

2, 7 (-^)'3/ J / + (2jF)> \09 + (')' O’ + e) s 4 t/ 7 (— /3) 5 / 4 

105 [( - (~/?) 1/2 ] 21 J, 15 .,[(-(-/ 9 ) 1 ^ 13 J; 

2r ) ' , (-fl) 1 / 1 V(-/3)”/ 4 + 2 tj 7 (-/3) 7 /' 

3»2[{ - M) 1 ' 2 ] _ 192., fajf - (-m i 

V(-|8) ,3/4 4i,(-/3 )»/< 2, : '(-,3) 1 ' 1 / 4 / ’ 

dr} 


V(~/5 ) 17 / 4 

(3-42) 




77 = ± 


e - (-/3) 1/2 1 


2*7 




« + f/3+|(-« 3/2 ]’ /, i f?(- l 9 ) ,/2 


(3.43) 


(3.44) 


It can be easily shown using the large argument form of the Fresnel transition 
function (see Appendix C) that when £ » \~py 12 or 77 > 1, Equation (3.37) 
appropriately reduces to the first six terms of Equation (3.4). Also, Equations (3.41) 
and (3.42) remain finite near the caustic when ( — ♦ (— /3 ) 1 / 2 and 77 — » 0 , however, 
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tern 


& 


they become numerically unstable and an alternative formulation should be used. 
Applying L’Hospital’s rule in Equation (3.18), and then using Equations (3.2) and 
(3. 19)— (3.32) we have: 


E{0,t,V * 0) ~ 


/o 


m-0) 


i 




+ 


b f 


9 o 


H 3 T /4 (- / 3 ) 3/2 

9 iV 


(2 j) 2 (-0Y /2 

hu 


1 - 


(-P) 31 ' + (-0f n J 

h\V h 2 T} 2 


929 


+ 


1- 


(-/3)3/4 T (-^9)3/2 


(2 m-Ty 

ko 

(2j) 4 (-/?) 11/z I" (~P) 3/4 "" (-/?) 3/2 J 


k iV 


+ 


W 


(3.45) 


and 


«o) *5^ 


-W l 


7/m 


Zl + 

4(-/3)V< 2 2(-/J)V» 


5/^ 


<7o 


*[£ - 


( 2 i ) 2 (-^) 7/2 

ft-0 


j\ 13firi 77 g x | 4g 2 T? 2 

j 1 4(-/?) 3 /< 2 + (-/3) 3 / 2_r (-0) 1 / 2 

f 19fcn/ 

1 4(- 


I 


(-pyn 


vm-er 

fcu f 25kiJi k, 
(~l n l 4( 3) J " 2 


(2j) 


_h 1 iim ! 

2 + 2(-/3 ) 3 / 2 

M£ - (-ff) 1/z ] 

(-/3) 3/! (-W /2 


) 


)• 


(3.46) 


where 

fo = 1/3, /, = 1.25, U = 4/3, 


</(j = 0.29629630, <71 = 3.00781250, g 2 — 5.83333333, 
hu = 0.69135802, h, = 4.74804688, h 2 = 13.3333333, 
Jfe 0 = 2.63374486, fc, = 6.48405151, fc 2 = 23.8333333. 


Also, A(t/) and its derivative near the caustic are given by: 


A (77 « 0) 


and 


^A(, « 0) 


2 




p jV 


-f-dW’l 1 / 2 





(3.47) 

(3.48) 
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Chapter 4 

Numerical Results and Discussion 

For an efficient and accurate computation of the incomplete Airy function gtj(/3,£) 
and its derivative, the argument space is divided into three regions as shown in 
Figure 4.1. Three different sets of formulae are used, one for each region in Figure 4.1, 
and a fourth set of formulae that is used in the immediate vicinity of the caustic 
(P + e 0 on; « 0). In region I, Equations (3.4) and (3.5) are used, in region 
II, Equations (3.37)-(3.44) are used, and in region III, the series solution is used 
given by Equations (2.21) and (2.22). In the immediate vicinity of the caustic and 
specifically when tj < 0.1, Equations (3.37)-(3.40) and (3.44)-(3.48) are used. 

An empirical expression for the number of terms used in the series solution is 
given by: 

JV(/3,0 = 8|/3| + 4, when £ < 2.0 , (4.1) 

— 8|/3| + 4 + 3|/3|(£ — 2.0) , when £ ^ 2.0 , (4.2) 

and results in a truncation error of less than 10 -6 . 

Figure 4.2 shows the percent amplitude error of the asymptotic result relative 
to the series solution along the boundary between regions I and III. The results are 
plotted vs. the parameter /3, with £ = (12 — 2/3) 1 / 2 . The asymptotic result shows 
excellent agreement with the series solution, exhibiting a maximum error of 0.12%. 
Figure 4.3 shows the percent amplitude error of the asymptotic result relative to the 
series solution along the boundary between regions II and III. The results for this 
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Figure 4.1: Three different sets of fo rm ulae are used for the computation of the 
incomplete Airy function, one for each region in the figure, and a fourth set that is 
used in the immediate vicinity of the caustic. 




0.0 1 .0 2J0 3J0 4J0 5.0 6.0 


P 

Figure 4.2: Percent amplitude error of the asymptotic result for the incomplete Airy 
function (solid line) and its derivative (broken line) along the boundary between 
regions I and III. Results are plotted vs. the parameter /3 with £ = (12 — 2/3)'^ 2 . 

case are plotted vs. the parameter £, with jd = —4. Again the asymptotic result 
shows excellent agreement with the series solution, exhibiting a maximum error of 
only 0.075%. 

Figures 4.4 and 4.5 show plots of the incomplete Airy function g 0 (/?,£) and 
its derivative vs. the parameter /? for £ = —3 and 2, respectively. Figures 4.6 
and 4.7 show plots of the incomplete Airy function g 0 (/?)0 and its derivative vs. the 
parameter ( for (3 = — 5 and 0, respectively. The marks on the contours represent 
direct numerical integration data, and show excellent agreement with the results 
obtained using the formulae derived in this report. 


24 





Figure 4.3: Percent amplitude error of the asymptotic result for the incomplete Airy 
function (solid line) and its derivative (broken line) along the boundary between 
regions II and III. Results are plotted vs. the parameter ( with /? = —4.0. 
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Figure 4.4: Plots of the real part (solid line) and imaginary part (broken line) of the 
incomplete Airy function and its derivative vs. the parameter (3 with £ = —3.0. The 
marks on the contours represent direct numerical integration data. 
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Figure 4.5: Plots of the real part (solid line) and imaginary part (broken line) of the 
incomplete Airy function and its derivative vs. the parameter with £ = 2.0. The 
marks on the contours represent direct numerical integration data. 
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Figure 4.6: Plots of the real part (solid line) and imaginary part (broken line) of the 
incomplete Airy function and its derivative vs. the parameter ( with /? = —5.0. The 
marks on the contours represent direct numerical integration data. 
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Chapter 5 

Summary and Conclusions 


In this report, a convergent series solution form for the incomplete Airy functions 
has been derived, and asymptotic expansions involving several terms have been 
developed for use in large argument approximations. It has been demonstrated that 
the combination of the series solution with the asymptotic formulae results in an 
efficient and accurate means for computing the incomplete Airy functions for the 
entire range of their arguments, 

A necessary requirement for the applicability of uniform asymptotic solutions 
to practical engineering problems is the efficient and accurate computation of the 
canonical functions involved. The results of this report would allow the formula- 
tion of useful uniform asymptotic solutions for several high-frequency scattering and 
diffraction problems in which the incomplete Airy integrals serve as canonical func- 
tions for the description of high-frequency field behavior in the vicinity of composite 
shadow boundaries. Furthermore, the methods used in this report may provide use- 
ful insight to the computation of other multivariable canonical functions occurring 
in high-frequency scattering and diffraction theory. 

A FORTRAN code for the computation of the incomplete Airy functions based 
on the formulae derived in this report is available from the authors. A complete 
code listing appears in Appendix D. 
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Appendix A 

Uniform Asymptotic Evaluation 
of Certain Stationary Phase 
Integrals 


Let’s consider the uniform asymptotic evaluation of a stationary phase integral of 
the form 

/(a, b; k) = r f(s)e jk ^' b Us (A.l) 

Ja 

where the phase function ) possesses two stationary phase points 5 ^ 2 ( 6 ) sa t- 

isfying 2 ,fc) = 0 with no restrictions placed on their location relative to the 
integration endpoint a . The amplitude function f(s) is assumed to be a slowly vary- 
ing and analytic function of s. The integral in (A.l) may be transformed into a 
canonical form using the following transformation: 


— 


4>(s, b ) = r(z,/?) = a -f fiz + z 3 / 3 

(A.2) 


where 

. ■ 


* : 


q = f(0,/3) = (f>(s p ,b ) ; <j>"{s p ) = 0 

0 3 

(A.3) 



e - * w [r'Mi 

(A.4) 


or alternatively 


= T(z h „f3)+1(-I3f*, 
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a 


(A.5) 


0 = -{^(0,/?)-r(* u ,/3)]} !, \ 


(A.6) 


with Z\ 2 = ±(— /3) 1 / 2 . The proper branch for /3 depends on the sign of <f>'{s,,) and 
<f>'"(s p ). Thus, using (A. 2) the integral in (A.l) becomes: 


/(&,,/?;*) = / °° G(z)e ikr ^ dz 


(A.7) 


where 




I 


6 = (« - * P ) 

0(z) = /(.)£, and 

<X5 


dz 


= [<£'(“) - ^'(*p )] 1/2 


1* 


r(*p)J 


<T(*p) 


(A.8) 
(A.9) 
(A- 10) 


The proper branches for and jt depend on the sign of <f>"\a p ). Next, we employ 
the Chester et. al. expansion [10] for the amplitude function in (A.7), i.e., 


c(*)= m (z i +Dr+b m i(z i +0r] 


(A. 11) 


m= 0 


and since only the leading terms in the asymptotic expansion of (A.7) will be re- 
tained, Equation (A. 11) may be rewritten as follows: 


(7(z) — a-o + zby + (z 2 + 


(A. 12) 


where 


a u 

bo 

ff( z ) 


+ G(z 2 ) , 

1 

2zi 

OO 


[G{z\ — G(z 2 ] , 

5^[a m (2 2 +/9) m_1 + b m z(z 2 + p)™*' 1 ] , and 

2/? 1 / 2 


m = l 




= /(*u) 


*=*1,2 


*"(»., j)J 


(A.13) 

(A.14) 

(A.15) 
(A. 16) 
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Inserting (A. 12 ) into (A. 7) yields 


1 


a 0 f“°°" e M0*+ z */ 3 ) dz -f 6 0 / °° ze ik V* +z * /3) dz 
J^a ''ta 


gtf a y k w*Hlw _ JL r°°” ^( z ) e ifc(/3*+* s /3) dz 

1 K J f a 


(A.17) 


jkJb 

The last two terms in the right side of (A.17) are obtained by an integration by parts 
of the integrand involving 5 ( 2 ). The desired uniform asymptotic approximation for 
It — ► 00 is given by the first three terms of (A.17), and may be expressed in the form 

I(U,m ~ e»« 

1 [ G(t,a) ~ OP ~ 61 fy) 


^fc(^Ca+^/3) 


) 


M e a 2 + /5 

where Ai(er,C) is the incomplete Airy integral defined by 


(A.18) 


li(cr,C)i r 'e^’/Vdz 


(A. 19) 


and the upper limit terminates in one of the three sectors in the complex z-plane 
where the integral converges. For example, in the case when the upper limit termi- 
nates in sector 0 < argz < 7t/3 we have that Aj(fr,£) = g 0 ( cr » 0- 


33 


Appendix B 

Computation of ao(£), &i(£) and 
their Derivatives 


The functions a 0 (0 and bi(() needed in the series solution form of the incomplete 
Airy functions can be expressed in terms of the incomplete Gamma function; i.e., 


a 0 (0 = e ^/6 3-2/3 r(1/3j _^3/3) j (B.l) 

6 t (0 = -e' W6 3 " 1/3 T(2/3, -yf/3), (B.2) 

where 

/»oc 

T(x,y) = t x - ] e- l dt, %{x)>0. (B.3) 

Jy 

Using the series solution form of the incomplete Gamma function [7], au(£) and i>i(£) 
are computed using the expressions: 

»«({) ~ ^ 6 3- 2/3 r(l/3)-£E(^7i^i. ( B ' 4 > 

and M£) - - e ' W 63 ' 1 , 3 r( 2 / 3 )-if 2 E ■ ( B5 > 

The number of terms in the series, N(£), is given by the empirical formula: 

N(£) = 2 + 4£ 2 , (B.6) 

and results in a truncation error of less than 10 -6 . 
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The derivatives of a ( i(£) and b t (£) are given by: 


•f’K) = 

jd k a { 0 k 3 \() + e fc aj,* 2) (0 + !) (0 , n > 3 , 

(B.T) 

II 

- l)af -,) (() 1 + j(a { 0 k) (() , n > 1 , 

(B.8) 

where 

ai'to = -e* V3 . 

(B.9) 


and ni 2, (f) = . 

(B.10) 

The constants and 

e k are obtained using the recursive relationships: 



d k = flfjfc-i + e*_i , 

(B.ll) 


and e k = e*_i + 2 , k >3 , 

(B.12) 


with c /2 = e 2 = 0. 
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Appendix C 

Computation of the Fresnel 
Transition Function 


The Fresnel transition function [3] is defined as follows: 

F(x) = 2 jyfi e jx e~ jr7 dr , 

where y/x = |-y/x| if * > 0, and y/x = j|-y/x| if x < 0. Also, 

When x < 6.0, F(x) is computed using its series solution form; i.e., 

, . Ar ( x ) (-j X ) n 

F(x) ~ y/irx e^ x+n ^' > sign(x) — 2jxe* x ^ 


(C.l) 


(C.2) 


(C.3) 


^ (2n + l)n! 

The number of terms in the series, N(x) } is given by the empirical formula: 

N(x) = 10y/i, (C.4) 

and results in a truncation error of less than 10 -6 . The large argument form of the 
Fresnel transition function (x > 6.0) is given by: 


. (-l) m 2 m (1), 
F(x) ~ 2^ 


771 — 0 




(C.5) 


where 


+ =1, and 2 m (a + ^) = (2a + l)(2a + 3) • • • (2a + 2k - 1) . (C.6) 
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Appendix D 
Code Listing 


C 


C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 
C! ! 


OPTIONS/EXTEND_SOURCE 

SUBROUTINE INCAIRY (B , X , GO , GOP , Gi , GiP , G2 , G2P) 

IMPLICIT NONE 

This subroutine computes the incomplete Airy functions and their 
first derivatives with respect to the parameter B. Both B and X 
are considered real. Double precision arithmetic is used. 


GO (B , X) = INT (X , INFTY ,FS) 

GOP (B ,X) = INT (X , INFTY , CJ*S*FS) 


where FS = EXP [CJ* (B*S+S**3/3)] 


Gl (B ,X) = GO(B,X)-TPI*AI(B) 

G2(B,X) = GO (B , X) -PI* [AI (B)+CJ*Bl(B)] 

AI (B) and BI(B) are the complete Airy functions. 


Version i.l 4-12-1992 

Author: E. D. Constantinides 

The Ohio State University 
ElectroScience Laboratory 
1320 Kinnear Road 
Columbus, OH 43212 

REAL *8 B , X , XP , X2 , BXX , PI/3 . 141592653589793/ , 

ft TPI/6 . 283185307179586/ 
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COMPLEX* 16 GO ,GOP ,G1 , G1P ,G2 ,G2P ,CJ/ (O.ODO, 1 .ODO)/ 
COMPLEX BZ,AI,AIP,BI,BIP 

EXTERNAL INC_AIRY_SS,INC_ATRY_FF,INC_AIRY_IP,AIBI 
COMMON CJ ,PI ,TPI 

C! ! ! 

XP=DABS (X) 

X2=X*X 

BXX=B+X2 

IF (B.GT. O.ODO) BXX=2 . ODO*B+X2 
IF (BXX.GT.12.ODO) THEN 

CALL INC_AIRY_IP (B , XP , GO , GOP) 

ELSEIF (B.LT.-4.0D0) THEN 

CALL INC_AIRY_FF (B , XP , GO , GOP) 

ELSE 

CALL INC_AIRY_SS (B ,XP ,GO ,GOP) 

ENDIF 

BZ=CMPLX(B) 

CALL AIBI (BZ , AI , AIP ,BI ,BIP) 

IF (X.LT. O.ODO) THEN 
GO=TPI*AI-DCONJG(GO) 

GOP=TPI*AIP-DCONJG(GOP) 

ENDIF 

Gl=GO-TPI*AI 

G1P=G0P-TPI*AIP 

G2=GO-PI*(AI+CJ*BI) 

G2P=G0P-PI* (AIP+C J*BIP) 

C! ! ! 

END 

C 

OPTIONS/EXTEND.SOURCE 

SUBROUTINE INC_AIRY_IP (BS , XS , GO , GOP) 

IMPLICIT NONE 

C! ! ! 

C! ! ! This subroutine computes the incomplete airy function 
C ! ! ! using integration by parts when BS+X5**2»1.0 
C! * ! 

REAL *8 BS ,XS ,BXS ,BXSi , BXS2 , BXS3 , BXS4.BXS5 ,BXS6 , BXS7 . 

ft BXS8,BXS9,BXS10,BXSli ,PI,TPI 

COMPLEX* 16 CJ,CT,GO,GOP 
COMMON CJ ,PI ,TPI 

cm 

BXS=BS+XS*XS 
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BXS1=1 . ODO/BXS 

BXS2=BXS1*BXS1 

BXS3=BXS2*BXS1 

BXS4=BXS3*BXS1 

BXS5=BXS4*BXS1 

BXS6=BXS5*BXS1 

BXS7=BXS6*BXS1 

BXS8=BXS7*BXS1 

BXS9=BXS8*BXSi 

BXS10=BXS9*BXS1 

BXS11=BXS10*BXS1 

CT=CDEXP(CJ* (BS*XS+XS*XS*XS/3 . ODO) ) 


IF (XS.EQ.O.ODO) GOTO 10 

GO=CJ*BXS1+2.0DO*XS*BXS3+2^0DO*CJ*BXS4-12.QDO*CJ*XS*XS*BXS5 

ft +40 . 0D0*XS*BXS6+40 . 0D0*CJ*BXS7-120 . 0D0*XS*XS*XS*BXS7 
ft -840 . 0D0*CJ*XS*XS*BXS8+1680 . ODO*C J*XS*XS*XS*XS*BXS9 

ft +2240 . 0D0*XS*BXS9 

GO=CT*GO 

GOP=-C J+BXS2-6 . 0D0*XS*BXS4-g . 0D0*CJ*BXS5+60 . 0D0*CJ*XS*XS+BXS6 
ft -240 . 0D0*XS*BXS7-280 .0D0*CJ*BXS8+840 . 0D0*XS*XS*XS*BXS8 

ft +6720.0D0*CJ*XS*XS*BXS9-15120.0D0*CJ*XS*XS*XS*XS*BXS10 

ft -20160. 0D0*XS*BXS10 

GOP=CJ*XS*GO+CT*GOP 


RETURN 

CM! = 

10 GO=C J*BXS 1+2 . 0D0*CJ*BXS4+4O . 0D0+CJ+BXS7 

GOP=-C J+BXS2-8 . ODO*C J+BXS5-280 . 0D0*CJ*BXS8-22400 . ODO*CJ*BXSU 
RETURN 

C! ! ! 

END 

C — 

OPTIONS/EXTEND.SOURCE 

SUBROUTINE INC_AIRY_FF(BS,XS,GO,GOP) 

IMPLICIT NONE 

C! ! ! 

C ! ! ! This subroutine computes the incomplete airy function 
C! ! ! using the Fresnel integral when BS<<0.0 
C! ! ! 

REAL *8 BS,XS,BSP,Z, ZZ ,BXS ,BX ,PI ,TPI ,B1 ,B2 ,B3 , 

ft G1 ,G2 , G3 ,Z3,Z5 ,Z7 ,Z9 , BXS1 ,BXS2 ,BXS3 , BXS4 , BXS5 ,BXS6 , 

ft BXS7 ,BXS8 ,BS14,BXi ,BSQP ,BSS ,F1 ,F2,F3 ,F4 ,F11 ,F12 ,F21 , 

ft F22 ,F31 ,F32 ,F41 ,F42 ,DZXS,DZBS 
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COMPLEX* 1 6 C J , CT , GO , GOP , FCTZ , CZZ , GOl , G02 , G03 , 
ft EPC , EPCP , SPC , SPCP , CT1 , GZ , FZ , HZ , KZ , FFCT , GPZ , FPZ , HPZ , 

ft KPZ.FCTZP 

EXTERNAL FTRANSD 
COMMON CJ ,PI ,TPI 


DATA 

ft 

DATA 

DATA 

DATA 

ft 

DATA 

ft 


G1,G2,G3 

FI , Fll ,F12 
F2,F21,F22 
F3 ,F3i ,F32 

F4,F41,F42 


/O . 2083333333333333 , 0 . 3342013888888889 , 
1.025812596450617/ 

/O . 3333333333333333 , 1 . 25D0 , 1 . 333333333333333/ 

/O . 2962962962962963 , 3 . 0078125D0 . 5 . 833333333333333/ 
/O . 6913580246913580 , 4 . 748046875D0 , 
13.33333333333333/ 

/2. 633744855967078, 6. 484051513671875, 

23 . 83333333333333/ 


BSP=DABS (BS) 

BSQP=DSQRT (BSP) 

BS14=DSqRT(BSqP) 

BSS=BSP*BSqP 

Bl=l . ODO/BSS 

B2=B1/BSS 

B3=B2/BSS 

BX=BS*XS+XS*XS*XS/3 . ODO 
BXl=-2 . 0D0+BSS/3 . ODO 
ZZ=DABS (BX+2 . 0D0*BSS/3 . ODO) 
Z=DSqRT(ZZ) 

IF (XS.LT.BSqP) Z=-Z 
Z3=Z*ZZ 


Z5=Z3*ZZ 


Z7=Z5*ZZ 


Z9=Z7*ZZ 

CT=CDEXP(CJ*BX) 

CT1=CDEXP(CJ*BX1) 

G01=0 . 5D0*CJ*G1*B1 
G02=-0 . 25D0*G2*B2 
G03=-0 . 125D0*CJ*G3*B3 
SPC=1 . 0D0+G01+G02+G03 

SPCP= ( 1 . 5D0*G01+3 . 0D0*G02+4 . 5D0*G03) /BSP 

BXS=BS+XS*XS 

DZXS=XS-BSqP 

CZZ=CDEXP(CJ*ZZ) 

IF (ABS(BXS) .LT.0.1D0) GOTO 10 
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iiUllili I i, J , hi* III I lilt U 11 


C! ! ! 


5 BXS1=1 . ODO/BXS 

— BXS2=BXS1*BXS1 

BXS3=BXS2*BXS1 
BXS4=BXS3*BXS 1 

- BXS5=BXS4*BXS1 

BXS6=BXS5*BXS1 

| BXS7=BXS6*BXS1 

BXS8=BXS7*BXS1 
CALL FTRANSD (ZZ , FFCT) 

FCTZ=0 . 5D0*CJ*C0NJG (FFCT) * CZZ/ Z 
FCTZP=-0 . 5DO*CZZ*DZXS/Z 
IF (Z.LT.O.ODO) FCTZ=FCT Z+CDSQ RT (CJ*PI) 

GZ=2 . ODO*BXSi-l . 0D0/Z/BS14 ~ 

GPZ=-2 . 0D0*BXS2-0 . 25D0/Z/BSP/BS14+0 . 5D0*DZXS/Z3/BS14 
FZ=-8 . 0D0*XS*BXS3+1 . 0D0/Z3/ BS 14-G1*B1/Z/BS14 
FPZ=24 . 0D0*XS*BXS4+0 . 25D0/Z3/BS14/BSP-1 .75D0*G1*B1/Z/BS14/BSP 
ft -1 . 5D0*DZXS/Z5/BS14+0 . 5D0*G1*B1*DZXS/Z3/BS14 

HZ=-16 .0D0*BXS4+96 .0D0*XS*XS*BXS5-3.0D0/Z5/BS14+G1*B1/Z3/BS14 
ft -G2*B2/Z/BS14 

HPZ=64 . 0D0+BXS5-480 . 0D0*XS*XS*BXS6-0 . 75D0/Z5/BS14/BSP+ 
ft 1 . 75D0*Gl*Bl/Z3/BS14/BSP-3 . 25D0*G2*B2/Z/BS14/BSP+ 

ft 7 . 5DO*DZXS/Z7/BS14-1.5DO*DZXS*G1*B1/Z5/BS14+ 

ft 0 . 5D0*DZXS*G2*B2/Z3/BS14 

KZ=640 . 0D0*XS*BXS6-1920 . 0DQ*XS*XS*XS*BXS7+15 . 0D0/Z7/BS14- 
ft 3 . 0D0*G1*B1/Z5/BS14+G2*B27Z3/BS14-G3*B3/Z/BS14 
KPZ=-3840 . 0D0*XS*BXS7+13440 . 0D0*XS*XS*XS*BXS8+3 . 75D0/Z7/BS14/BSP 
ft -5 . 25D0*Gl*Bl/Z5/BS14/BSP+3. 25D0*G2*B2/Z3/BS 14/BSP 

ft -4 . 75*G3*B3/Z/BS14/BSP-52 . 5D0*DZXS/Z9/BS14 

ft +7 . 5D0*DZXS*Gl*Bl/Z7/BSi4-l . 5D0*DZXS*G2*B2/Z5/BS14 

ft +0.5D0*DZXS*G3*B3/Z3/BS14"" 

EPC=CT* (0 . 5D0*CJ*GZ-0 . 25D0*FZ-0 . 125D0*CJ*HZ+0 .0625D0*KZ) 

EPCP=CT* (0 . 5D0*CJ*GPZ-0 . 25D0+FPZ-0 . 125D0+CJ*HPZ+0 . 0625D0*KPZ) 
G0=CT1*SPC*FCTZ/BS14+EPC 

GOP= (CJ*BS14+0 . 25D0/BS14/BSP) *CT1*SPC*FCTZ+CT1*SPCP*FCTZ/BS14 
ft +CT1 *SPC*FCTZP/BS 14+C J *XS*EPC+EPCP 

RETURN 

C! ! ! 

CM! Small argument Form 
C! ! ! 

10 FCTZ=0.5D0*CDSQRT(CJ*PI)-Z-CJ*ZZ*Z/3.0D0 

DZBS=0 . 5D0/BS 14/DSQRT ( 1 . ODO+DZXS/3 . ODO/BSQP) 
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FCTZP=-DZBS*CZZ 

GZ=1 . 0D0-F11*Z/BSQP/BS14+F12*ZZ/BSS 

GPZ=-0.75D0*Fll*Z/BSqP/BS14/BSP+1.5D0*F12*ZZ/BSS/BSP 
ft -Fll*DZBS/BSqP/BS14+2 . 0D0*Z*F12*DZBS/BSS 

FZ=1 .0D0-F21*Z/BSQP/BS14+F22*ZZ/BSS 

FPZ=- 0 .75D0*F21*Z/BSqP/BS14/BSP+i . 5D0*F22*ZZ/BSS/BSP 
ft -F21*DZBS/BSqP/BS14+2 .0D0*Z*F22*DZBS/BSS 

HZ=1 . 0D0-F3i*Z/BSqP/BS14+F32*ZZ/BSS 

HPZ=-0 . 75D0*F31*Z/BSqP/BS14/BSP+l . 5D0*F32*ZZ/BSS/BSP 
ft -F31*DZBS/BSqP/BS14+2 . 0D0*Z*F32*DZBS/BSS 

KZ=1 . 0D0-F41*Z/BSqP/BS14+F42*ZZ/BSS 

KPZ=-0 . 75D0*F41*Z/BSqP/BS14/BSP+i . 5D0*F42*ZZ/BSS/BSP 
ft -F41*DZBS/BSqP/BS14+2 . 0D0*Z*F42*DZBS/BSS 

EPC=CT* (-0 . 5*CJ*Fl*GZ/BSP+0 . 25D0*F2*FZ/BSP/BSS 
ft +0 . 125D0*F3*CJ*HZ/BSP/BSS/BSS 

ft -0 . 0625D0*F4*KZ/BSP/BSS/BSS/BSS) 

EPCP=CT* (-0 . 5*CJ*Fl*GZ/BSP/BSP-0 . 5*CJ*F1*GPZ/BSP 
ft +2 . 5D0*0 . 25D0*F2*FZ/BSP/BSS/BSP+0 . 26D0*F2*FPZ/BSP/BSS 

ft +4 . 0D0*0 . 125D0*F3*CJ*HZ/BSP/BSS/BSS/BSP 

ft +0 . 125D0*F3*CJ*HPZ/BSP/BSS/BSS 

ft -5 . 5D0*0 . 0625D0*F4*KZ/BSP/BSS/BSS/BSS/BSP 

ft -0 . 0625D0*F4*KPZ/BSP/BSS/BSS/BSS) 

G0=CT1*SPC*FCTZ/BS14+EPC 

G0P= (C J*BS 14+0 . 25D0/BS14/BSP) *CT1*SPC*FCTZ 
ft +CTl*SPCP*FCTZ/BS14+CTi*SPC*FCTZP/BS14+Cj*XS*EPC+EPCP 

C ! ! ! 

RETURN 

END 

C 

OPTIONS/EXTEND.SOURCE 
SUBROUTINE FTRANSD (XF ,FFCT) 

IMPLICIT NONE 

cm 

C! ! ! This routine evaluates the Fresnel Transition function F(X) 

C! ! ! 

REAL *8 X ,XF ,PI4 ,PI ,TPI 

INTEGER N.NT 

COMPLEX* 16 C J , AO , AN , FFCT , FFCTS 
COMMON CJ,PI,TPI 

C! ! ! 

X=DABS (XF) 

IF (X.GT.6.0D0) GOTO 1 
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C! ! ! 

C! ! ! Small argument (series) form 

cm ' 

IF (X.Eq.O.ODO) THEN 
FFCT= (0 . 0D0 ,0 .0D0) 

RETURN 
END IF 

PI4=PI/4.0D0 

FFCT=DSqRT(PI*X)*CDEXP(CJ*(X+PI4)) 

A0= ( 1 . 0D0 , 0 . 0D0) 

FFCTS=A0 

NT=10*DSqRT(X) 

DO N=1 ,NT 

AN=-CJ*X*AO/DBLE(N) 

FFCTS=FFCTS+AN/DBLE (2*N+iJ 
A0=AN 
END DO 

FFCT=FFCT-2 . 0D0*C J*X*CDEXP (CJ*X) *FFCTS 
GOTO 20 

C! ! ! 

C! ! ! Large argument form 
CM! 

1 A0=(1 . 0D0.0.0D0) 

FFCT=A0 
DO N=1 ,8 

AN=0 . 5DO*CJ*DBLE(2*N-l)*AO/X 
FFCT=FFCT+AN 
A0=AN 
END DO 

20 IF (XF.GE.O.O) RETURN 

FFCT=DC0N JG (FFCT) 

C! ! ! 

RETURN 

END 

C 

OPTIONS/EXTEND_SOURCE 

SUBROUTINE INC_AIRY_SS(BS,XS,GO,GOP) 

IMPLICIT NONE 

C! ! ! 

C ! !! This subroutine computes the Incomplete Airy function using a 
C! ! ! convergent series solution with error less than 1.0E - 6. 

C! ! ! 
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REAL*8 BS ,BSP ,BS2 ,BS3 , XS , XS2 ,XS3 ,C0 ,C1 ,PI ,TPI 

INTEGER M,N ,MM,NT,MT 

C0MPLEX*16 GO, GOP, A (66, 33) ,B(66,33) ,A0P(33) ,A1P(33) , 
ft CA(66) ,CB(66) ,CJ 

EXTERNAL AOA1XS 

COMMON CJ ,PI ,TPI 

C! ! ! 

CM! Initialize the coefficient values. 

C! ! ! 

BSP=DABS (BS) 

XS2=XS*XS 

XS3=XS2*XS 

CALL A0A1XS(XS ,AOP(l) ,A1P(1)) 

A0P(2)=-CDEXP (CJ*XS3/3 . ODO) 

AOP (3)=CJ*XS2*A0P(2) 

AOP (4) =2 . ODO*C J *XS*AOP ( 2 ) +C J *XS2*AOP (3) 
A1P(2)=CJ*XS*A0P(2) 

DO N=3 ,4 

A1P (N)=CJ*DBLE(N-2) *AOP(N-l)+CJ*XS*AOP (N) 

END DO 

BS2=BS*BS 

BS3=BS2*BS 

GO=AOP ( 1 ) +A1P ( 1 ) *BS+0 . 5DO*CJ*AOP (2 ) *BS2 
G0P=A1P(1) +CJ*AOP (2) *BS 
IF (BS.EQ.O.ODO) RETURN 
MT=8*BSP+4 

IF (XS . GT . 2 . ODO) MT=MT+3*BSP* (XS-2) 

NT=(MT+l)/2 
C0=0 . ODO 
Cl=2 . ODO 
DO M=1 ,3 
DO N=1 ,NT 

A(M,N)=(0. ODO , 0 . ODO ) 

B(M,N)=(0.0D0,0. ODO ) 

END DO 
END DO 

C! ! ! 

A(l,l)=(i. ODO , 0 . ODO) 

B(2,l)=(i. ODO ,0 .ODO) *BS 
A (3 , 2)=(0 .ODO ,0 . 5D0) *BS2 

C! ! ! 

C ! ! ! Compute the series 
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C! ! ! 



y 


i=i 

1 s 


DO M=4 , MT 
MM=(M+l)/2 

A (M , 1 ) =A (M-3 , 1) *BS3/DBLE(M-l) /DBLE(M-2) 

B (M , 1 ) =B (H-3 , 1 ) *BS3/DBLE (M- 1 ) /DBLE (M-2 ) 

DO N=2 ,MM 

A (M , N) = (A (M-3 , N) *BS3+C J*A (M-2 , N-l ) *BS2) /DBLE (M- 1 ) /DBLE (M-2 ) 
B(M,N)=(B(M-3 ,N)*BS3+CJ*B (M-2, N-l)*BS2) /DBLE (M-l) /DBLE (M-2) 
END DO 

DO N=MM+1 ,NT 

A (M , N) = (0 . ODO , 0 . ODO) 

B(M,N)=(O.ODO,O.ODO) 

END DO ~ , 

IF ((MM.NE. (M/2)) .AND. (MMTGE 5)) THEN 
C0=C0+C1 
Cl=Cl+2 .ODO 

A0P(MM)=CJ*XS2*A0P(MM-iy+Cl*CJ*XS*A0P(MM-2)+C0*CJ*A0P(MM-3) 
A1P (MM) =CJ*DBLE (MM-2) *AOP (MM-1 ) +C J*XS*AOP (MM) 

END IF 

CA (M) = (0 . ODO , 0 . ODO) 

CB (M) = (0 . ODO , 0 . ODO) 

DO N=i ,MM 

CA (M) =CA (M) +A (M , N) * AOP (N) 

CB (M) =CB (M) +B (M , N) * A1P (N) 

END DO 

GO=GO+(CA(M)+CB(M) ) 

GOP=GOP+DBLE (M- 1 ) * (CA (M) +CB (M) ) /BS 
END DO 

C! ! ! 

RETURN 

END 

- 

C 

OPTIONS/EXTEND.SOURCE 
SUBROUTINE A0A1XS(XS,A0,A1) 

IMPLICIT NONE 


REAL *8 XS , G13/2 . 6789385347D0/ ,G23/1 . 3541179394D0/ , PI ,TPI , 

ft XS 1 , XS2 , AIO/O . 3550280539D0/ , AIP0/-0 . 2588194038D0/ 

COMPLEX+16 A0,A1,CJ,X3,B(100),C13,C23 
INTEGER I, NT 

COMMON CJ.PI.TPI 
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DATA C13,C23/(0. 4163415888278022, 0.2403749283845681) , 

ft (-0 . 6004684775880014,0. 3466806371753174) / 

C! ! ! 

XS1=DABS (XS) 

X3=CJ*XS1*XS1*XS1/3.0D0 

XS2=XS1*XS1 

C! ! ! 

C! ! ! Series Solution 
C! ! ! 

A0=C13*G13 

Ai=C23*G23 

IF (XSl.Eq.O.ODO) RETURN 
B(1)=X3 

AO=AO-XS1-0. 25D0*XS1*B(1) 

Al=Al-0 . 5D0*CJ*XS2-0 . 2*C J*XS2*B (1) 

NT=2+4*XS2 
DO 1=2, NT 

B(I)=B(I-1) *X3/DBLE(I) 

A0=A0-XS1*B(I)/DBLE(3*I+1) 

Al=Al-CJ*XS2*B(I)/DBLE(3*I+2) 

END DO 

C! ! ! 

10 IF (XS.LT.O.ODO) THEN 

A0=2 . ODO*PI*AIO-DCONJG (AO) 

Al=2 . 0D0*PI*AIP0-DC0NJG(A1) 

END IF 

C! ! ! 

RETURN 

END 

C 

SUBROUTINE AIBI(Z,AI ,AIP,BI,BIP) 

C! ! ! 

C ! ! ! This routine calculates the Airy functions AI(Z),BI(Z), 

C ! ! ! and their derivatives AIP(Z) ,BIP(Z) . 

CM! Ref. Abramowitz and Stegun, Handbook of Mathematical Functions. 
C ! ! ! For CABS(Z) .LE. 6.0 ,a Taylor Series is used. 

C! ! ! ARG(Z) may take any value. See eqs . (10.4.2) to (10.4.5) . 

C! ! ! 

COMPLEX Z , AI , AIP ,BI ,BIP 

C! ! ! 

IF(CABS(Z) .GT.6.)G0 TO 12 
CALL AIBI1 (Z , AI ,AIP,BI ,BIP) 
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RETURN 

12 CALL AIBI2(Z,AI.AIP,BI,BI P ) 

C! ! ! 

RETURN 

END 

C 

SUBROUTINE AIBI1(Z,AI,AIP,BI,BIP) 

C! ! ! 

COMPLEX Z,AI,AIP,BI,BIP 
C0MPLEX*16 F,6,FP,GP 
DOUBLE PRECISION CC1.CC2 

C! ! ! 

DATA S3, CC1.CC2/1. 732050808, .355028053887817, .258819403792807 / 

C! ! ! 

CALL FZGZ(Z,F,G,FP,GP) 

AI=CCi*F-CC2*G 
AIP=CC1*FP-CC2*GP 
BI=S3* (CCi*F+CC2*G) 

BIP=S3* (CC1*FP+CC2*GP) 

C! ! ! 

RETURN 

END 

C 

SUBROUTINE FZGZ(Z ,F ,G ,FP ,GP) 

cm 

C ! f ! The auxiliary functions F(Z) ,G(Z) ,FP(Z) ,GP(Z) are computed as in 
C! ! ! "Tables of the Modified Hankel functions of order one-third and 
C ! ! ! of their derivatives". Computation Lab, Harvard Univ. PresB,1945. 
CM ! 

COMPLEX* 16 F , G , FP , GP , Z3 , Z3M, ZD 


COMPLEX 

Z 

REAL *8 

AM, BM, CM, DM, AO, BO, CO, DO 

REAL 

ZMBD(5) 

INTEGER 

MAX (5) 

DATA 

ZMBD /6.1, 5.6, 4.8, 4.1, 

DATA 

MAX /22, 19, 16, 14, 11 / 


ZD=0.D0 
ZD=Z 
A0=1 . DO 
B0=1 .DO 



1= 

LI 

U- 



C! ! ! 


C0=0 . 5D0 
D0=1 . DO 

Z3=(ZD**3)/200 
Z3M=Z3 
ZMAG=CABS (Z) 

DO 3 M=1 , 5 

3 IF(ZMAG .LE. ZMBD(M) )MADHAX=MAX(M) 

F=AO 

G=BO 

FP=CO 

GP=DO 

DO 10 M=1 ,MADMAX 
TM=FL0AT(3*M) 

AM=200 . DO*AO/TM/ (TM-1 ) 

BH=200 . DO*BO/TM/ (TM+1 ) 

CM=200 . DO*CO/TM/ (TM+2) 

DM=200 . DO*DO/TM/ (TM-2) 

F=F+AM*Z3M 

G=G+BM*Z3M 

FP=FP+CM*Z3M 

GP=GP+DM*Z3M 

Z3M=Z3M*Z3 

AO=AH 

BO=BM 

CO=CM 

DO=DM 

10 CONTINUE 
G=ZD*G 

FP=(ZD**2) *FP 

C! ! ! 

RETURN 

END 

C 

SUBROUTINE AIBI2 (XX , AI , AIP ,BI ,BIP) 

C! ! ! 

C! ! ! This Routine calculates the Airy functions AI(XX) ,BI(XX) , 

C! ! ! and their derivatives AIP (XX) ,BIP(XX) . 

C ! ! * Ref. Abramowitz and Stegun, Handbook of Mathematical Functions. 
C! ! ! 

COMPLEX Z, AI , AIP ,BI ,BIP ,XX 
COMPLEX Z25 , ZTB , ZT , ZT2 , ZT3 , ZT4 , ZT5 
COMPLEX CT1 , A2L2 , EIPI3 , EIPI6 , C , S 
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C! ! ! 


DATA RTPI ,TWORPI,RTOP,POF 

'/. /1. 772453851, 3. 544907702,. 797884561,. 785398164 / 

DATA A2L2.EIPI6.EIPI3 

’/, / (0 . , . 346573590) , ( . 8 6 6025404 , . 5) , ( . 5 , . 866025404) / 

DATA Cl/ . 069444444/ , C2/. 037133487/ ,C3/ . 037993059/ , 

'/, C4/ . 057649 190/, C5/ . 116099064/ 

DATA Dl/- . 097222222/ , D2/-. 043885030/ ,D3/- . 042462830/ , 

D4/~ . 062662163/ ,D5/- . 124105896/ 

C! ! ! 

ZTB=(2./3.)*XX**1.5 
ARG=ATAN2 (AIMAG (XX) , REAL (XX) ) 

IF (ABS (ARG) . GE . 2 . 1 ) GO TO 100 

CM! ~T,.' ■ 

CM! EQN. (10.4.59) , (10.4.61) 

C! ! ! 

Z25=XX** . 25 

ZT=ZTB 

ZT2=ZT*ZT 

ZT3=ZT2*ZT 

ZT4=ZT2*ZT2 

ZT5=ZT3*ZT2 

CT1=CEXP ( -ZT) /TWORPI 

AI=CT1/Z25*(1-C1/ZT+C2/ZT2-C3/ZT3+C4/ZT4-C5/ZT5) 

AIP=-CT1*Z25*(1-D1/ZT+D2/ZT2-D3/ZT3+D4/ZT4-D5/ZT5) 

IF (ARG .LT.O. )G0 TO 20 
ZT=(0. ,-l.)>fZTB 

C! ! ! 

CM! EQN. (10. 4. 65), (10. 4. 68) WITH UPPER SIGNS. 

C! ! ! 

Z=XX/EIPI3 
CT1=ZT+P0F-A2L2 
BI=EIPI6 
BIP=1 . /EIPI6 
GO TO 30 

20 ZT=(0 . , 1 . )*ZTB 

C! ! ! 

C! ! ! EQN. (10.4.65) , (10.4.68) WITH LOWER SIGNS. 

C! ! ! 

Z=XX*EIPI3 
CT1=ZT+P0F+A2L2 
BI=1 . /EIPI6 


49 



BIP=EIPI6 
30 S=CSIN (CT1) 

C=CC0S (CT1) 

Z25=Z**.25 

ZT2=ZT*ZT 

ZT3=ZT2*ZT 

ZT4=ZT2*ZT2 

ZT5=ZT3*ZT2 

BI=BI*RT0P/Z25*(S*(1-C2/ZT2+C4/ZT4)-C*(C1/ZT-C3/ZT3+C5/ZT5)) 
BIP=BIP*RT0P*Z25* (C* ( 1-D2/ZT2+D4/ZT4) +S* (D1/ZT-D3/ZT3+D5/ZT5) ) 
RETURN 

100 ZT=(0 . , 1 . )*ZTB 

C! ! ! 

CM! EQN. (10. 4. 60), (10. 4. 62), (10. 4. 64), (10. 4. 67) 

C! ! ! 

IF (ARC . LT . 0 . ) ZT=-ZT 
Z=-XX 

Z25=Z**.25 

ZT2=ZT*ZT 

ZT3=ZT2*ZT 

ZT4-ZT2*ZT2 

ZT5=ZT3*ZT2 

CT1=ZT+P0F 

S=CSIN(CT1) 

C=CC0S (CT1) 

AI=1 . /RTPI/Z25* (S* (l-C2/ZT2+C4/ZT4)-C*(Cl/ZT-C3/ZT3+C5/ZT5) ) 
AIP=-Z25/RTPI* (C* (1-D2/ZT2+D4/ZT4) +S* (D1/ZT-D3/ZT3+D5/ZT5) ) 
BI=1 . /RTPI/Z25* (C* (1-C2/ZT2+C4/ZT4)+S*(C1/ZT-C3/ZT3+C5/ZT5) ) 
BIP=Z25/RTPI* (S*(l-D2/ZT2+D4/ZT4)-C* (D1/ZT-D3/ZT3+D5/ZT5) ) 

C! ! ! 

RETURN 

END 

C 
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